###Propensity Score Simulation (slide 41-43)###

library(survey)

###Creat expit function###
expit<-function(x){
return(
exp(x)/(1+exp(x))
)
}

####Generate data###
set.seed(1000)
n<-5000
Z<-rbinom(n,1,.5)
C<-rnorm(n,Z,1)
A<-rbinom(n,1,expit(C-.5))
Y<-2*exp(C)+A+rnorm(n,0,1)

###Unadjusted difference in means###
summary(lm(Y~A))

###Adjusted difference in means (misspecified model)###
summary(lm(Y~A+C))

###Compute propensity scores###
p<-glm(A~C,family=binomial)$fitted.values

###Compute propensity scores weights###
wt<-A/p+(1-A)/(1-p)

###Create survey design object (no clusters, weights=wt)###
mydesign <- svydesign(ids = ~1, data = data.frame(Y,A), weights = ~wt)

###Propensity score weighted difference in means###
summary(svyglm(Y ~ A, design = mydesign))

###Adjusted difference in means (correctly specified model)###
expC<-exp(C)
summary(lm(Y~A+expC))


###Propensity Score Example (slide 45-60)###

library(mdscore)
library(multcomp)
library(survey)

###Load data###
booster<-read.table("C:\\Penndata\\class\\681\\booster.dat")
 
inj<-booster[,1] 
injnoc<-booster[,2] 
boost<-booster[,3] 
age_yrs<-booster[,4] 
vehicle<-booster[,5] 
youngd<-booster[,6] 
towaway<-booster[,7] 
intru<-booster[,8] 
airbag<-booster[,9]
car<-as.numeric(vehicle==0)
lvan<-as.numeric(vehicle==1)
truck<-as.numeric(vehicle==2)
mvan<-as.numeric(vehicle==3)
suv<-as.numeric(vehicle==4)

###t-tests for mean difference###
t.test(age_yrs[boost==0],age_yrs[boost==1])
t.test(car[boost==0],car[boost==1])
t.test(lvan[boost==0],lvan[boost==1])
t.test(truck[boost==0],truck[boost==1])
t.test(mvan[boost==0],mvan[boost==1])
t.test(suv[boost==0],suv[boost==1])
t.test(youngd[boost==0],youngd[boost==1])
t.test(towaway[boost==0],towaway[boost==1])
t.test(intru[boost==0],intru[boost==1])
t.test(airbag[boost==0],airbag[boost==1])

###Fit initial propensity score model###
summary(glm(boost~as.factor(age_yrs)+lvan+truck+mvan+suv+youngd+towaway+intru+airbag,family=("binomial")))

###Consider 2-way interactions###
base<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag,family=("binomial"))

ageveh<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
as.factor(age_yrs)*truck+as.factor(age_yrs)*mvan+as.factor(age_yrs)*suv,family=("binomial"))

lr.test(base,ageveh)

ageyoung<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
as.factor(age_yrs)*youngd,family=("binomial"))

lr.test(base,ageyoung)

agetow<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
as.factor(age_yrs)*towaway,family=("binomial"))

lr.test(base,agetow)

ageair<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
as.factor(age_yrs)*airbag,family=("binomial"))

lr.test(base,ageair)

vehyoung<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
youngd*truck+youngd*mvan+youngd*suv,family=("binomial"))

lr.test(base,vehyoung)

vehtow<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
towaway*truck+towaway*mvan+towaway*suv,family=("binomial"))

lr.test(base,vehtow)

vehair<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
airbag*truck+airbag*mvan+airbag*suv,family=("binomial"))

lr.test(base,vehair)

youngtow<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
youngd*towaway,family=("binomial"))

lr.test(base,youngtow)

youngair<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
youngd*airbag,family=("binomial"))

lr.test(base,youngair)

towair<-glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
towaway*airbag,family=("binomial"))

lr.test(base,towair)


###Update propensity score model to include age*vehicle and age*airbag interactions###
summary(glm(boost~as.factor(age_yrs)+truck+mvan+suv+youngd+towaway+airbag+
as.factor(age_yrs)*truck+as.factor(age_yrs)*mvan+as.factor(age_yrs)*suv+as.factor(age_yrs)*airbag,family=("binomial")))

###Rerun dropping age 7 * airbag interaction###
age5<-as.numeric(age_yrs==5)
age6<-as.numeric(age_yrs==6)
age7<-as.numeric(age_yrs==7)

summary(glm(boost~age5+age6+age7+truck+mvan+suv+youngd+towaway+airbag+
age5*truck+age5*mvan+age5*suv+
age6*truck+age6*mvan+age6*suv+
age7*truck+age7*mvan+age7*suv+
age5*airbag+age6*airbag,family=("binomial")))

###Final propensity score model###
propmod<-glm(boost~age5+age6+age7+truck+mvan+suv+youngd+towaway+airbag+
age5*truck+age5*mvan+age5*suv+
age6*truck+age6*mvan+age6*suv+
age7*truck+age7*mvan+age7*suv+
age5*airbag+age6*airbag,family=("binomial"))

propscore<-propmod$fitted.values

inj1<-inj[boost==1]
inj0<-inj[boost==0]

###Means and SEs of injury rates by booster seat status###
mean(inj1)
sd(inj1)/sqrt(length(inj1))
mean(inj0)
sd(inj0)/sqrt(length(inj0))

###Unadjusted difference in injury rates###
summary(lm(inj~boost))

###Weighted propensity score analysis###
###Create a vector of inverse propensity weights###
w <- boost / propscore + (1-boost) / (1-propscore)

###Create a survey design object###
mydesign <- svydesign(ids = ~1, data = data.frame(inj, boost), weights = ~w)

###Fit a weighted regression model###
mymodel <- svyglm(inj ~ boost, design = mydesign)

# Print summary of the weighted model
summary(mymodel)

###Stratified propensity score analysis###
###Find quintiles of propensity score###
cut<-sort(propscore)[8000*c(.2,.4,.6,.8)]

q1<-as.numeric(propscore<cut[1])
q2<-as.numeric((propscore>=cut[1]&(propscore<cut[2])))
q3<-as.numeric((propscore>=cut[2]&(propscore<cut[3])))
q4<-as.numeric((propscore>=cut[3]&(propscore<cut[4])))
q5<-as.numeric(propscore>=cut[4])

###T-test for first quintile of PS###
t.test(age_yrs[(boost==0)&(q1==1)],age_yrs[(boost==1)&(q1==1)])
t.test(car[(boost==0)&(q1==1)],car[(boost==1)&(q1==1)])
t.test(lvan[(boost==0)&(q1==1)],lvan[(boost==1)&(q1==1)])
t.test(truck[(boost==0)&(q1==1)],truck[(boost==1)&(q1==1)])
t.test(mvan[(boost==0)&(q1==1)],mvan[(boost==1)&(q1==1)])
t.test(suv[(boost==0)&(q1==1)],suv[(boost==1)&(q1==1)])
t.test(youngd[(boost==0)&(q1==1)],youngd[(boost==1)&(q1==1)])
t.test(towaway[(boost==0)&(q1==1)],towaway[(boost==1)&(q1==1)])
t.test(intru[(boost==0)&(q1==1)],intru[(boost==1)&(q1==1)])
t.test(airbag[(boost==0)&(q1==1)],airbag[(boost==1)&(q1==1)])

###T-test for second quintile of PS###
t.test(age_yrs[(boost==0)&(q2==1)],age_yrs[(boost==1)&(q2==1)])
t.test(car[(boost==0)&(q2==1)],car[(boost==1)&(q2==1)])
t.test(lvan[(boost==0)&(q2==1)],lvan[(boost==1)&(q2==1)])
t.test(truck[(boost==0)&(q2==1)],truck[(boost==1)&(q2==1)])
t.test(mvan[(boost==0)&(q2==1)],mvan[(boost==1)&(q2==1)])
t.test(suv[(boost==0)&(q2==1)],suv[(boost==1)&(q2==1)])
t.test(youngd[(boost==0)&(q2==1)],youngd[(boost==1)&(q2==1)])
t.test(towaway[(boost==0)&(q2==1)],towaway[(boost==1)&(q2==1)])
t.test(intru[(boost==0)&(q2==1)],intru[(boost==1)&(q2==1)])
t.test(airbag[(boost==0)&(q2==1)],airbag[(boost==1)&(q2==1)])

###T-test for third quintile of PS###
t.test(age_yrs[(boost==0)&(q3==1)],age_yrs[(boost==1)&(q3==1)])
t.test(car[(boost==0)&(q3==1)],car[(boost==1)&(q3==1)])
t.test(lvan[(boost==0)&(q3==1)],lvan[(boost==1)&(q3==1)])
t.test(truck[(boost==0)&(q3==1)],truck[(boost==1)&(q3==1)])
t.test(mvan[(boost==0)&(q3==1)],mvan[(boost==1)&(q3==1)])
t.test(suv[(boost==0)&(q3==1)],suv[(boost==1)&(q3==1)])
t.test(youngd[(boost==0)&(q3==1)],youngd[(boost==1)&(q3==1)])
t.test(towaway[(boost==0)&(q3==1)],towaway[(boost==1)&(q3==1)])
t.test(intru[(boost==0)&(q3==1)],intru[(boost==1)&(q3==1)])
t.test(airbag[(boost==0)&(q3==1)],airbag[(boost==1)&(q3==1)])

###T-test for fourth quintile of PS###
t.test(age_yrs[(boost==0)&(q4==1)],age_yrs[(boost==1)&(q4==1)])
t.test(car[(boost==0)&(q4==1)],car[(boost==1)&(q4==1)])
t.test(lvan[(boost==0)&(q4==1)],lvan[(boost==1)&(q4==1)])
t.test(truck[(boost==0)&(q4==1)],truck[(boost==1)&(q4==1)])
t.test(mvan[(boost==0)&(q4==1)],mvan[(boost==1)&(q4==1)])
t.test(suv[(boost==0)&(q4==1)],suv[(boost==1)&(q4==1)])
t.test(youngd[(boost==0)&(q4==1)],youngd[(boost==1)&(q4==1)])
t.test(towaway[(boost==0)&(q4==1)],towaway[(boost==1)&(q4==1)])
t.test(intru[(boost==0)&(q4==1)],intru[(boost==1)&(q4==1)])
t.test(airbag[(boost==0)&(q4==1)],airbag[(boost==1)&(q4==1)])

###T-test for fifth quintile of PS###
t.test(age_yrs[(boost==0)&(q5==1)],age_yrs[(boost==1)&(q5==1)])
t.test(car[(boost==0)&(q5==1)],car[(boost==1)&(q5==1)])
t.test(lvan[(boost==0)&(q5==1)],lvan[(boost==1)&(q5==1)])
t.test(truck[(boost==0)&(q5==1)],truck[(boost==1)&(q5==1)])
t.test(mvan[(boost==0)&(q5==1)],mvan[(boost==1)&(q5==1)])
t.test(suv[(boost==0)&(q5==1)],suv[(boost==1)&(q5==1)])
t.test(youngd[(boost==0)&(q5==1)],youngd[(boost==1)&(q5==1)])
t.test(towaway[(boost==0)&(q5==1)],towaway[(boost==1)&(q5==1)])
t.test(intru[(boost==0)&(q5==1)],intru[(boost==1)&(q5==1)])
t.test(airbag[(boost==0)&(q5==1)],airbag[(boost==1)&(q5==1)])

fit<-lm(inj~q2+q3+q4+q5+boost+boost*q2+boost*q3+boost*q4+boost*q5)
summary(fit)

###Compute ACE from stratified model###
###Compute fraction of sample in each propensity score quintile###
N<-8000
P<-c(sum(q1==1),sum(q2==1),sum(q3==1),sum(q4==1),sum(q5==1))/N
summary(glht(fit,linfct=t(c(0,0,0,0,0,1,P[2],P[3],P[4],P[5]))))

###Propensity Score Overlap Example (slide 71-76)###

library(survey)

booster<-read.table("C:\\Penndata\\class\\681\\booster.dat")
 
inj<-booster[,1] 
injnoc<-booster[,2] 
boost<-booster[,3] 
age_yrs<-booster[,4] 
vehicle<-booster[,5] 
youngd<-booster[,6] 
towaway<-booster[,7] 
intru<-booster[,8] 
airbag<-booster[,9]
car<-as.numeric(vehicle==0)
lvan<-as.numeric(vehicle==1)
truck<-as.numeric(vehicle==2)
mvan<-as.numeric(vehicle==3)
suv<-as.numeric(vehicle==4)

###Normalized distance###

nd<-function(vari)
{
m1<-mean(vari[boost==1])
m0<-mean(vari[boost==0])
s21<-var(vari[boost==1])
s20<-var(vari[boost==0])
delta<-(m1-m0)/sqrt((s21+s20)/2)
return(delta)
}

nd(age_yrs)
nd(car)
nd(lvan)
nd(truck)
nd(mvan)
nd(suv)
nd(youngd)
nd(towaway)
nd(intru)
nd(airbag)

###Overall measure###
tX<-cbind(age_yrs,lvan,truck,mvan,suv,youngd,towaway,intru,airbag)
tX1<-tX[(boost==1),]
tX0<-tX[(boost==0),]
meandiff<-apply(tX1,2,mean)-apply(tX0,2,mean)
var<-(cov(tX1)+cov(tX0))/2
sqrt(t(meandiff)%*%solve(var)%*%meandiff)

###Generate propensity score again###
age5<-as.numeric(age_yrs==5)
age6<-as.numeric(age_yrs==6)
age7<-as.numeric(age_yrs==7)

propmod<-glm(boost~age5+age6+age7+truck+mvan+suv+youngd+towaway+airbag+
age5*truck+age5*mvan+age5*suv+
age6*truck+age6*mvan+age6*suv+
age7*truck+age7*mvan+age7*suv+
age5*airbag+age6*airbag,family=("binomial"))

propscore<-propmod$fitted.values

###Box plot of PS by exposure###
boxplot(propscore~boost)

###Range of PS by exposure###
min(propscore[boost==0])
max(propscore[boost==0])

min(propscore[boost==1])
max(propscore[boost==1])

###N of children not overlapping###
sum(propscore<.0248)

###Summary of child factors for those not overlapping###
drop<-(propscore<.0248)
summary(age_yrs[drop])
table(vehicle[drop])/94
table(youngd[drop])/94
table(towaway[drop])/94
table(airbag[drop])/94

###Redo PS-weighted analysis dropping these children###
###Set weight to 0 for non-overlapping children###
nw<-w*(propscore>=.0248)

mydesign<-svydesign(ids = ~1, data = data.frame(inj,boost), weights = ~nw)
mymodel<-svyglm(inj~boost,design=mydesign)
summary(mymodel)


###Matching example (slides 90-101)###

###Malhalanobis distance function###
malhal<-function(x0,x1){
return(
sqrt(
(x0-x1)%*%Siginv%*%(x0-x1)
)
)
}

booster<-read.table("C:\\Penndata\\class\\681\\booster.dat")
 
inj<-booster[,1] 
injnoc<-booster[,2] 
boost<-booster[,3] 
age_yrs<-booster[,4] 
vehicle<-booster[,5] 
youngd<-booster[,6] 
towaway<-booster[,7] 
intru<-booster[,8] 
airbag<-booster[,9]
car<-as.numeric(vehicle==0)
lvan<-as.numeric(vehicle==1)
truck<-as.numeric(vehicle==2)
mvan<-as.numeric(vehicle==3)
suv<-as.numeric(vehicle==4)

age5<-as.numeric(age_yrs==5)
age6<-as.numeric(age_yrs==6)
age7<-as.numeric(age_yrs==7)

###Compute propensity score###
propmod<-glm(boost~age5+age6+age7+truck+mvan+suv+youngd+towaway+airbag+
age5*truck+age5*mvan+age5*suv+
age6*truck+age6*mvan+age6*suv+
age7*truck+age7*mvan+age7*suv+
age5*airbag+age6*airbag,family=("binomial"))

###Sort data by propensity score###
propscore<-propmod$fitted.values

###Create temporary design matrix to sort by PS
tX<-cbind(age_yrs,lvan,truck,mvan,suv,youngd,towaway,intru,airbag)

###Sort by PS###
X<-tX[order(propscore,decreasing=TRUE),]
Y<-inj[order(propscore,decreasing=TRUE)]
BOOSTER<-boost[order(propscore,decreasing=TRUE)]

###Create design matrix and outcome stratified by exposure###
X1<-X[(BOOSTER==1),]
X0<-X[(BOOSTER==0),]
Y1<-Y[(BOOSTER==1)]
Y0<-Y[(BOOSTER==0)]

###Compute MH distance###
mu1<-apply(X1,2,mean)
mu0<-apply(X0,2,mean)
Sig1<-cov(X1)
Sig0<-cov(X0)
Sig<-.5*Sig1+.5*Sig0
Siginv<-solve(Sig)

N<-8000
N1<-sum(boost)
N0<-N-N1

holdj<-rep(0,N1)
holdtest<-rep(0,N1)
X0<-cbind(c(1:N0),X0)
nX0<-X0

###Generate distance measures for each observation###
for(i in 1:N1){
test<-rep(0,(N0-(i-1)))
for(j in 1:(N0-(i-1))){
test[j]<-malhal(nX0[j,2:length(nX0[1,])],X1[i,])
}
###Find control (seat belt) with closest MH distance to treated (booster seat)###
holdtest[i]<-min(test)
holdj[i]<-nX0[order(test)[1],1]
j<-order(test)[1]
nX0<-nX0[-j,]
print(i)
}

plot.ts(holdtest,main="MH Distance by Order of PS")

###Compute matched ACE for treated cases###
tau<-mean(Y1-Y0[holdj])

vtau<-sum((Y1-Y0[holdj]-tau)^2)/(N1*(N1-1))

tau
c(tau-1.96*sqrt(vtau),tau+1.96*sqrt(vtau))

###Check balance###
mean(X1[,1])-mean(X0[,2])
mean(X1[,1]-X0[holdj,2])

mean(X1[,7])-mean(X0[,8])
mean(X1[,7]-X0[holdj,8])

mean(X1[,9])-mean(X0[,10])
mean(X1[,9]-X0[holdj,10])

###Compute ACE for treated and control cases and combine###
###Compute tau-hat-t###
diff1<-rep(0,N1)
M1<-rep(0,N1)
for(i in 1:N1){
test<-rep(0,N0)
for(j in 1:N0){
test[j]<-malhal(X0[j,2:length(X0[1,])],X1[i,])
}
M1[i]<-sum(test<=1)
diff1[i]<-Y1[i]-sum(Y0[test<=1])/sum(test<=1)
print(i)
}

tau1<-mean(diff1[M1>0])
s12<-var(Y1)
s02<-var(Y0)
vtau1<-(s12+s02/mean(M1[M1>0]))/sum(M1>0)

###Compute tau-hat-ctau1###
diff0<-rep(0,N0)
M0<-rep(0,N0)
for(i in 1:N0){
test<-rep(0,N1)
for(j in 1:N1){
test[j]<-malhal(X0[i,2:length(X0[1,])],X1[j,])
}
M0[i]<-sum(test<=1)
diff0[i]<-sum(Y1[test<=1])/sum(test<=1)-Y0[i]
print(i)
}

tau0<-mean(diff0[M0>0])
s12<-var(Y1)
s02<-var(Y0)
vtau0<-(s02+s12/mean(M0[M0>0]))/sum(M0>0)

###Obtain full ACE estimator###
tau<-(N0/(N0+N1))*tau0+(N1/(N0+N1))*tau1
vtau<-((N0/(N0+N1))^2)*vtau0+((N1/(N0+N1))^2)*vtau1
c((tau-1.96*sqrt(vtau)),(tau+1.96*sqrt(vtau)))

###Repeat using boostrap###
BOOT<-100
bootau<-rep(0,BOOT)
for(boot in 1:BOOT){

resamp1<-sample(c(1:N1),replace=TRUE)
resamp0<-sample(c(1:N0),replace=TRUE)

tX0<-X0[resamp0,]
tX1<-X1[resamp1,]
tY0<-Y0[resamp0]
tY1<-Y1[resamp1]

diff1<-rep(0,N1)
M1<-rep(0,N1)
for(i in 1:N1){
test<-rep(0,N0)
for(j in 1:N0){
test[j]<-malhal(tX0[j,2:length(tX0[1,])],tX1[i,])
}
M1[i]<-sum(test<=1)
diff1[i]<-tY1[i]-sum(tY0[test<=1])/sum(test<=1)
#print(i)
}

tau1<-mean(diff1[M1>0])

diff0<-rep(0,N0)
M0<-rep(0,N0)
for(i in 1:N0){
test<-rep(0,N1)
for(j in 1:N1){
test[j]<-malhal(tX0[i,2:length(tX0[1,])],tX1[j,])
}
M0[i]<-sum(test<=1)
diff0[i]<-sum(tY1[test<=1])/sum(test<=1)-tY0[i]
#print(i)
}

tau0<-mean(diff0[M0>0])

bootau[boot]<-(N0/(N0+N1))*tau0+(N1/(N0+N1))*tau1

print(boot)
}

var(bootau)

###Repeat for propensity score###

propdist<-function(p1,p2){
return(
(log(p1/(1-p1))-log(p2/(1-p2)))^2
)
}

propscoresort<-propmod$fitted.values[order(propscore,decreasing=TRUE)]

propscoresort1<-propscoresort[(BOOSTER==1)]
propscoresort0<-propscoresort[(BOOSTER==0)]


###Check balance###
holdj<-rep(0,N1)
holdtest<-rep(0,N1)
id0<-c(1:N0)

for(i in 1:N1){
test<-rep(0,(N0-(i-1)))
for(j in 1:(N0-(i-1))){
test[j]<-propdist(propscoresort0[j],propscoresort1[i])
}
holdtest[i]<-min(test)
holdj[i]<-id0[order(test)[1]]
print(i)
}

mean(X1[,1])-mean(X0[,2])
mean(X1[,1]-X0[holdj,2])

mean(X1[,7])-mean(X0[,8])
mean(X1[,7]-X0[holdj,8])

mean(X1[,9])-mean(X0[,10])
mean(X1[,9]-X0[holdj,10])

mean(X1[,1])-mean(X0[,2])
mean(X1[,1]-X0[holdj,2])

mean(X1[,7])-mean(X0[,8])
mean(X1[,7]-X0[holdj,8])

mean(X1[,9])-mean(X0[,10])
mean(X1[,9]-X0[holdj,10])

###Compute ACE for treated and control cases and combine###
###Compute tau-hat-t###
diff1<-rep(0,N1)
M1<-rep(0,N1)
for(i in 1:N1){
test<-rep(0,N0)
for(j in 1:N0){
test[j]<-propdist(propscoresort0[j],propscoresort1[i])
}
M1[i]<-sum(test<=.05)
diff1[i]<-Y1[i]-sum(Y0[test<=.05])/sum(test<=.05)
print(i)
}

tau1<-mean(diff1[M1>0])
s12<-var(Y1)
s02<-var(Y0)
vtau1<-(s12+s02/mean(M1[M1>0]))/sum(M1>0)

###Compute tau-hat-c###
diff0<-rep(0,N0)
M0<-rep(0,N0)
for(i in 1:N0){
test<-rep(0,N1)
for(j in 1:N1){
test[j]<-propdist(propscoresort0[i],propscoresort1[j])
}
M0[i]<-sum(test<=.05)
diff0[i]<-sum(Y1[test<=.05])/sum(test<=.05)-Y0[i]
print(i)
}


tau0<-mean(diff0[M0>0])
s12<-var(Y1)
s02<-var(Y0)
vtau0<-(s02+s12/mean(M0[M0>0]))/sum(M0>0)

tau<-(N0/(N0+N1))*tau0+(N1/(N0+N1))*tau1
vtau<-((N0/(N0+N1))^2)*vtau0+((N1/(N0+N1))^2)*vtau1
c((tau-1.96*sqrt(vtau)),(tau+1.96*sqrt(vtau)))

###Bootstrap###
BOOT<-100
bootau<-rep(0,BOOT)
for(boot in 1:BOOT){

resamp1<-sample(c(1:N1),replace=TRUE)
resamp0<-sample(c(1:N0),replace=TRUE)

tpropscoresort0<-propscoresort0[resamp0]
tpropscoresort1<-propscoresort1[resamp1]
tY0<-Y0[resamp0]
tY1<-Y1[resamp1]

diff1<-rep(0,N1)
M1<-rep(0,N1)
for(i in 1:N1){
test<-rep(0,N0)
for(j in 1:N0){
test[j]<-propdist(tpropscoresort0[j],tpropscoresort1[i])
}
M1[i]<-sum(test<=.05)
diff1[i]<-Y1[i]-sum(Y0[test<=.05])/sum(test<=.05)
#print(i)
}

tau1<-mean(diff1[M1>0])

diff0<-rep(0,N0)
M0<-rep(0,N0)
for(i in 1:N0){
test<-rep(0,N1)
for(j in 1:N1){
test[j]<-propdist(tpropscoresort0[i],tpropscoresort1[j])
}
M0[i]<-sum(test<=.05)
diff0[i]<-sum(Y1[test<=.05])/sum(test<=.05)-Y0[i]
#print(i)
}

tau0<-mean(diff0[M0>0])

bootau[boot]<-(N0/(N0+N1))*tau0+(N1/(N0+N1))*tau1

print(boot)
}

var(bootau)
